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Abstract. Previously, an explicit solution for the time evolution of the Wigner 
function was presented in terms of auxiliary phase space coordinates which obey simple 
equations that are analogous with, but not identical to, the classical equations of 
motion. They can be solved easily and their solutions can be utilized to construct the 
time evolution of the Wigner function. In this paper, the usefulness of this explicit 
solution is demonstrated by solving a numerical example in which the Wigner function 
has strong spatial and temporal variations as well as regions with negative values. It 
is found that the explicit solution gives a correct description of the time evolution 
of the Wigner function. We examine next the pseudoparticle method which uses 
classical trajectories to evolve the Wigner function. We find that the lowest-order 
pseudoparticle approximation reproduces the general features of the time evolution, 
but there are deviations. We show how these deviations can be systematically reduced 
by including higher-order correction terms in powers of Ti 2 . 



f Invited talk presented at the Wigner Centennial Conference, Pecs, Hungary, July 8-12, 2002, to be 
published in the Journal of Optics B: Quantum and Classical Optics, June 2003. 
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1. Introduction 

It is a great pleasure for me to participate in the Wigner Centennial Conference at 
Pecs, Hungary in honor of the Centennial of Professor Eugene Wigner's birthday. I 
first met Professor Wigner in 1959 when I was an undergraduate student at Princeton. 
Professor Wigner taught us and our fellow graduate school classmates Curtis Callen, 
Stephen Adler, and Alfred Goldhaber a course in Advanced Quantum Mechanics in 
1961. We were at the celebration party at the Graduate College when Professor Wigner 
was awarded the Nobel Prize in Physics in 1963. After I obtained my Ph. D. degree 
at Princeton, I came to work at Oak Ridge National Laboratory, for which Professor 
Wigner was the founding Scientific Director. I saw Professor Wigner from time to 
time at Oak Ridge when he came to work on his civil defense project. It is therefore 
particularly meaningful to me to come and celebrate his centennial in his native land to 
commemorate his many important contributions. 

As a part of the "Junior Paper" research project at Princeton University, I went 
to see Professor Wigner in his office in Fine Hall in 1959 and inquired about the early 
history of Quantum Mechanics. He told me that the great thing about Schrodinger was 
that his wave equation was formulated in configuration space. However, to Professor 
Wigner, who was trained as a chemical engineer, dynamics could also be described 
in phase space. His formulation of quantum mechanics in phase space in 1932 led 
to the well-known Wigner description of quantum systems pQ. The joint function 
of coordinate and momentum f(rp) introduced in 1932, now known as the Wigner 
function, is analogous with the classical distribution function. This analogy has provided 
new insights and useful applications to many quantum systems. It is deservedly a subject 
of intense interest in this Wigner Centennial Conference [2]. 

The Wigner function is, however, not identical to the classical distribution function. 
The classical distribution function is always non-negative and can be described in terms 
of a collection of cell points with positive weights in phase space. The evolution of 
the classical distribution function can be followed by tracking these cell points using 
classical equations of motion [3]. 

The Wigner function can assume negative as well as positive values in some regions 
of phase space. How does a Wigner function with regions of negative values evolve as 
a function of time? The concept of a probability distribution function with positive 
weights is clearly not applicable here. The difficulty of propagating the Wigner function 
with negative values has been a barrier to the study of the dynamics of the Wigner 
function using the equation of motion of the Wigner function directly. 

Previously, I formulated a simple method to propagate the Wigner function in time 

The ability to propagate the Wigner function with negative values arises from 
the fact that one deals with amplitudes in this method, and one does not need to use 
the concept of probabilities. The method contains three important components. First, 
one employs an auxiliary variable s, which is allowed to span the whole configuration 
space. Second, using this auxiliary variable s, the Wigner function can be represented 
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in terms of auxiliary phase space coordinates R and P. The equations of motion for 
these auxiliary coordinates -R and P are quite simple. They are analogous with, but 
not identical to, the classical equations of motion, as dP/dt depends on the auxiliary 
coordinate s. They can be solved easily and their solutions can be utilized to construct 
the Wigner function at the next time step. The third component consists of providing 
the correct amplitude function for the evolution. For each value of s and the evolution 
of phase space coordinates, one associates an amplitude factor exp{is ■ (p — P)/h}. 
After the amplitudes for all possible values of s have been obtained at the end of one 
time step, one adds all these amplitudes together coherently (Huygen Principle) to give 
the Wigner function f(rp) at the next time step. This explicit solution overcomes the 
difficulty of propagating a Wigner function with negative values. 

Recently, with the rapid advances in quantum information technology and the 
reduction in the size of micro-electronic devices, there has been renewed interest in 
experimental and theoretical studies of the Wigner function 4j- 20J. Experimental 
measurements of the Wigner function have been carried out using many different 
techniques Regions of negative Wigner function were observed in many 

experiments IB]- The time evolution of the Wigner function has been studied by 
many authors |TU] — [21]. It is instructive to review the explicit solution for the time 
evolution of the Wigner function obtained previously and to demonstrate its usefulness 
by a numerical example. We choose to examine an example in which the time evolution 
of the Wigner function can be ready evaluated by using wave functions and eigenvalues. 
The Wigner function in the example should have significant spatial and temporal 
variations as well as regions of negative values to test whether the explicit solution 
obtained previously is capable of treating these non-classical features. The success in 
providing a correct time evolution of the Wigner function using this explicit solution 
will pave the way for its applications in the future. 

As an approximation of the explicit solution, a pseudoparticle method was 
developed previously [3] to obtain the time evolution of the Wigner function by following 
the trajectories of the phase space coordinates, using classical equations of motion . It is 
instructive to review this pseudoparticle method and assess its accuracy by comparing 
with the correct results. Furthermore, it is useful to develop systematic ways to improve 
the accuracy of the pseudoparticle method to assist its future applications. 

2. Explicit Solution of the Time Evolution of the Wigner Function 

We shall first review the explicit solution obtained in 1982 3\. We consider a particle 
in the potential V 




(1) 



The Wigner function is 




(2) 
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From the Schrodinger equation, we obtain the equation of motion for the Wigner 
function [I] 

df{rp,t) p 



V P /(rp,f) - |sin ||v r y ■ V^j V(r,t)f(rp,t) = 0, 



(3) 



dt m 

where Vjf acts on the potential V and acts on the Wigner function /. 

Given the Wigner function f(r Q p Q ,t ) at time to, we wish to obtain the Wigner 
function at the next time step at t = to + dt with a small St. To solve for the dynamics, 
we divide the phase space into cells (pseudoparticles). Consider one such cell centered 
at {r p } at time to with volume element dr dp and amplitude f(rop ,to). We 
retain the auxiliary coordinate s and follow the time dependence of each cell in the 
Lagrangian sense. For the phase space coordinates {^oPol initially at time to, we 
label their phase space coordinates at subsequent time t by {R(rop s, t) P(roP Q s,t)}. 
In classical dynamics, the evolution of the phase space coordinates will not depend 
on s. In quantum dynamics, they will depend on s. The momentum coordinate P 
can jump to the momentum coordinate p at time t. We associate an amplitude of 
exp{is • [p — P(roP s, t)]/h} for this momentum jump. This amplitude factor is chosen 
because for the classical case where P is independent of s this amplitude leads to the 
correct delta function propagator 5(p — P(r p ,t)) after we integrate over s [see Eq. 
(JHJ)]. We therefore express the Wigner function f{rp,t) at the next time step as |3j 
' dr dp ds 
(27rft) J 

Our task is to find solutions for the phase space coordinates R(r p s, t) and P(r p s, t) 
with initial conditions -R(r p s,t ) = r and P(r p s,to) = Po- If we can nn d these 
coordinates at the next time step, the above integral can then be carried out to give the 
Wigner function at the new time t. 

To obtain the equations of motion for R(r p s,t) and P(r p s, t), we substitute 
Eq. (jlj) into Eq. (jSJ). For the first term in Eq. Q, we get 
df(rp,t) f dr dp ds 



-jyLexp{is- \p-P(r p s,t)}/h}5[r-R(r p s,t)]f(roP ,t ). (4) 



dt J (2vrft) 3 
-is) dP 



-f(roP ,t- 



exp{is ■ [p - P(r p s, t)\/%}5[r - R(r p s, t)] 



h dt 

+ exp{is ■ [p - P(r p s, t)]/h}(V R 5[r - R(r p s, t)]) ■ dR(r p s, t)/dt .(5) 
Noting that 

V R 6[r - R(ro Po s, t)} = -V r S[r - P(r p s, *)], (6) 

we can rewrite the second term inside the square bracket as 
•dr dp ds 

7^F /(roPo ' to) 

x exp{is ■ [p - P(r PoS,t)]/h}(-WrS[r - R(r p s,t)] ■ dR(r p s,t)/dt 
dR{r p s,t) 



dt 



V r f(rp,t), (7) 
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where we have used Eq. to obtain the right-hand side. For the last term in Eq. (jSJ), 
we substitute (HJ) into this term and find 

h. 




VJT • V p '} V(r,t)f(r P ,t) = U^0^f(roPo,to) 



\V(r - - t) - V(r + - t)] 
x exp{zs • [p- P(r oPo s,t)]/h}5[r - R(r oPo s,t)]) [ -^ 2 ' } - l ^ 2 ' JJ . (8) 

Putting Eqs. ©-© into Eq. fl3J), we obtain 

aR(w, t) +j ,i 

at mj 



rdr dp ri ds „. 



(27rft) 



is) «9P [^(r-f,*)-V(r + f,*)] 



x exp{is ■ \p - P(r p s, t)]/h}5[r - R(r p s, t)) = 0. (9) 

By comparing similar terms, the above equation leads to the equations of motion for R 
and P 0, 

f - (io) 

ot m 

In the above equations, r and R are interchangeable, because of the 5-function 5(r — R) 
in Eq. It is interesting to note that these two equations were used recently by 

Morawetz to study the quantum response of a finite fermion system |19j . 

Equation (jlj) combined with the equations of motion (fTUj) and (|TT|) gives an explicit 
and exact solution of the time evolution of the Wigner function [3]. The equation of 
motion for R is the same as the classical equation of motion. The equation of motion 
for P coincides with the classical equation of motion for small s but differs from the 
classical equation when s is large. The latter difference distinguishes a quantum system 
from a classical system. 

We note that if (1) we expand Eq. (fTT]) in s and we keep only the lowest order 
term, or if (2) we take the limit h — > so that the dominant contributions in Eq. 
come from regions of small s, then Eqs. (ITUJ) and (JTTJ) becomes the classical equations 
of motion [3| 

ot m 

and 

^ = -V R V(R,t). (13) 
The function P(r p s,t) is then independent of s, and Eq. (jH) gives 

f(rp, t) = J dr dp 5[r - R(r p , t)}5[p - P(r p , t)]f(r p , t ), (14) 
which is the same as solution of the evolution of the classical distribution function. 
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3. Implementation of the Explicit Solution 



The explicit solution of Eqs. (jlj), (fTUJ) and (fTTJ) given in the previous section allows one 
to propagate the Wigner function from one time to another. One considers a small time 
increment St = t — to, and solves Eqs. (fT0|) and |TT| . One obtains 



P 

r = R(r p s,t) = r H St. 

m 



and 



V r 



Po 



2 ' 



St, 



(15) 



(16) 



where e s is the unit vector in the s direction. After r(R) and P for different s have 
been obtained, we can substitute them in Eq. (jlj) and carry out the integration to give 
the new Wigner function at the next time step. Eqs. (fTKJ) and (fTKJ) lead to 



f(rp,t) = dr dp t 



ds 



(27Th) : 



exp 



St' 

J 



p 



xS(r -r St)f(r ,p ,t ). 



rn 



(17) 



We shall first consider the one-dimensional case for which the factor (27th) 3 becomes 
2irh. The integration over Xo(^o) can be easily carried out. the above equation for the 
one-dimensional case becomes 



f(xp,t) = dp 2 



00 ds 



o 2nh 
p 



cos 



s ■ (p - po) 



h 



V(x- S -,t)-V(x+ S -,t) 



St 

J 



*f(x St,p ,t ). 

m 



This explicit solution of the time evolution of the Wigner function can be written as 



f(xp,t) = J dp ^F c (x,p -po) + F s (x,p-p )\f(x - ^St,p ,to), 

where F c (x,p — po) and F s (x,p — po) are cosine and sine transforms given by 



F c (x,p -p ) = 2 
and 

F s (x,p -p ) = 2 



00 ds 



o 2nh 

00 ds 



o 2nh 



cos 



sin 



s ■ (p ~ Po) 

h 



s-(p~ Po) 

h 



COS 



sin 



V(x- S -,t)-V(x+ S -,t) 



V(x- S -,t)-V{x+ S -,t) 



St 

~h 
St 



(19) 



(20) 



(21) 



The above set of equations (fT9*|) - ()21|) provide a simple implementation of the explicit 

solution (Jlj), (fTTj). and (fTUj) of the time evolution of a one-dimensional Wigner function. 

If one keeps only terms up to first order in St for small values of St, the Wigner 

function evolution equation becomes 

p St r p 

f(xp,t)=f(x St,p,t ) + — dp V s (x,p-p )f(x St,p ,t ), (22) 

m n J m 

where 

V s (x,p-p ) = 



ds 



(2-nhy 



sin 



s-(p~ Po) 

h 



V(x- S -,t)-V(x+ S -,t)\. 



(23) 



Explicit Solution of the Time Evolution of the Wigner Function^ 7 

We shall now discuss the three-dimensional case for which the integration over 
s in Eq. (|T7|) is rather complicated. It can be simplified if one makes the following 
approximation 

V(r - |, t) - V(r + 1, f) « £ (V(r - S -f-,t)- V(r + t)) , (24) 

where the righthand side contains all terms of the lefthand side except the cross terms 
involving products of different Sj and Sj>(j' ^ j) with high-order derivatives of V(r), 
when one expand the above functions in powers of Sj. If one uses this approximation 
by neglecting these cross terms, Eq. (|T7|) for the three-dimensional case becomes 



f(rp,t)= d Po H 2 

J 3=1 



00 dsi 



o 2nh 



cos 



s j(Pj -Poj) 

h 



V(r 



st 
h 



p 



x/(r 5t,p ,t ). 

m 



(25) 



The explicit solution of the time evolution of the Wigner function can be written as 
f(rp,t) = j dp Y[(F cj (r, Pj - p 0j ) + F sj (r, Pj - p 0j )\f{r - —St,p Q ,t Q ), (26) 
where F C j(r,pj — poj) and F S j(r,pj — poj) are cosine and sine transforms given by 



F cj (r, Pj -p oj ) = 2 
and 

F sj{r,Pj -Poj) = 2 



00 dsj 



o 27ih 



00 dsi 



cos 



Sj{Pj -Poj) 



COS 



V( r -^ l ,t)-V(r + 



s j e j 



,t) 



5t 
J 
(27) 



o 2nh 



sin 



h 



sin 



V( r -^,t)-y(r + ^,t))£ 



St 

h 

(2; 



The above set of equations ()26 p -()28| l provide a simple and approximate implementation 
of the explicit solution (0J), (jllj) . and of the time evolution of a three-dimensional 
Wigner function. 



4. An Example of the Evolution of the Wigner Function 

It is of interest to consider an explicit example to test whether the solutions in Eqs. (|19jl- 
(|2*T|) or (J22J)-((2HJ) leads to the correct results. For this purpose we study the dynamics 
of the Wigner function for a case that can be easily evaluated and compared with the 
Wigner function obtained by using the explicit solution. 

We examine a particle in a one-dimensional attractive Gaussian potential well. We 
express physical quantities in dimensionless units, with a Schrodinger equation given by 

mx) = {-l^~ e ~ X2/2a2 } nx) - (29) 
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The wave function can be expanded in terms of a set of non-orthogonal Gaussian 
wave functions with different widths, 

*(ar) = ]T a n ip n (x), (30) 

n=l 

where the normalized basis function (x\n) = ip n (x) is taken to be 

e -x 2 /2f3l 

Mx) = (31) 

and /3 2 = u(3q. The Hamiltonian matrix in this basis set can be easily constructed. 
Specifically, one calculates the matrix element of the overlap matrix B for this set of 
basis states 



Bnm = (n\m) = )j-jj!rfjg< (32) 
the matrix element of the kinetic energy matrix T 



{nl \ w) |m> -W#+/M+/%' (33) 



| v 1 / WrA 

2dx 2 \ m 2V^ + / 
and the matrix element for potential energy matrix V 

V nm = (n\(-e- x2 ^ 2 )\m) 

where 



2P n (3 m a 2 

\ Q2 , Q2 Q2 — 1— a 2 



(34) 



B 2 B 2 

y-nm Q2 \ Q2 ' V ' 

The eigenvalue equation becomes 

(T + V)a = EBa, (36) 

where a is the column matrix of the coefficients {a n } of Eq. 1)30)1 . This eigenvalue 
equation can be diagonalized to yield the eigenenergies 1*^ and eigenfunctions ^a^) 
represented by the set of coefficients {a\ n }. We normalize ^\(x) by / dx\^\(x)\ 2 = 1. 

A simple non-stationary state $(a;) of the system can be constructed as a linear 
combination of the eigenstates ^\(x) with amplitudes b\, 

$(x,t = 0) = J2bxVx(x). (37) 

A 

The time dependence of this non-stationary system is 

$(x,t)=J2he- tE ^x(x). (38) 

A 

The Wigner function of the system is then given by 

f(xp,t) = Y, c n(ty m (t)f nm (xp), (39) 

nm 

where 

Cn(t) = E & Ae-^a A „, (40) 
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fnmM = 2,l^^e-^, (41) 



in 



and 



Vnm — m m ■, V^^) 



(43) 




Unm — 4 m t 2 X I ] + tp 



2$ 

In our numerical example, we study the dynamics of a system in a potential with a 
width parameter a = 3. The two lowest energy even-parity eigenstates ^oO^) an d *ffi(x) 
have eigenvalues E = —0.844 and E 1 = —0.312 respectively. As a test, we construct a 
non-stationary wave function with equal amplitudes in these two eigenstates at t — 0, 

$(x, t = 0) = 4|{*o(^) + (44) 

The time dependence of the wave function is then 

$( x , t ) = -^{e- iEot y (x) + e-^V^x)}. (45) 
v2 

The time dependence of the Wigner function can be evaluated using Eq. (j3Hjl . The 
Wigner functions for different values of p and t are shown in Figs, (la), (lb), and 
(lc) on the left panel of Fig. 1. The Wigner function has both positive and negative 
values in different regions of phase space. The eigenf unctions and eigenvalues give a 
correct representation of the Wigner function as it depends only on the accuracy of the 
eigenvalues and eigenf unctions, and not on the accuracy of the method of propagating 
forward in time. 

The non-stationary state we have used leads to an initial Wigner function with 
significant spatial and temporal variations as well as regions of negative values. These 
properties provide a stringent test of the methods of using Eqs. (fT§j) - (fTTj) or Eqs. 
(j2HJ) to propagate the Wigner function. 

For our tests we take the initial Wigner function for our non-stationary state Eq. 
f!44jl at time t = 0, and evolve the Wigner function using the explicit solution. We carry 
out the time evolution using small time increments. We can choose either Eqs. (|19j)- 
121(1 or Eqs. (|22 )) -(|23 )) to evolve the Wigner function. The former includes higher order 
effects in 5t, and can be used with a larger time increment than the second approach. We 
choose to use Eqs. (fTT?|) - (j2T]) . and the results for the time evolution from t — to t — 3 
in Figs, (la'), (lb'), and (lc') were calculated using 30 time steps. The cosine and sine 
transforms can be calculated using the method of fast Fourier transforms. The numerical 
integration over po in Eq. ()19)) can then be carried out to yield the Wigner function at 
the next time step. This process is repeated in a stepwise manner to propagate forward 
in time. The resulting time evolution of the Wigner function is shown in Figs, (la'), 
(lb'), and (lc') on the right panel of Fig. 1. 

A comparison of the left and the right panels indicates that the explicit solution 
of Eqs. fll9|) - (f2*T| gives an excellent reproduction of the time evolution of the Wigner 
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Figure 1. Time dependence of the Wigner function obtained in two different methods. 
The left panel gives the Wigner function obtained from eigenfunctions and eigenvalues, 
and the right panel shows that Wigner function obtained by the explicit solution of 
Eqs. P) -ft Hl . 



function obtained by using eigenvalues and eigenfunctions. The phases and the negative 
regions of the Wigner function are correctly reproduced. The difference between the 
eigenfunction solution and the explicit solution is small. For example, for p = 0.6 at 
t = 3, the maximum value of the eigenfunction solution of the Wigner function is 0.9313 
at x — 1.692, and the explicit solution gives 0.9206 at x — 1.692. For this momentum 
p and t = 3, the minimum of the eigenfunction solution of the Wigner function is 
—0.3888 at x = —1.1667, and the the explicit solution gives —0.3735 at x = —1.1667. 
The positions of the maxima and minima in the two methods are the same, and the 
magnitudes differ by about 1%. The explicit solution of Eqs. fl"9|) - (f2T|) thus leads to an 
accurate determination of the time evolution of the Wigner function, even if it contains 
regions of negative or oscillating values. 
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5. Pseudoparticle Method and Approximations 

As an approximation to obtain the explicit solution of Eqs. (J2J), (fTUJ) and (JTTJ), a 
pseudoparticle method was presented previously to calculate the time evolution of the 
Wigner function 0. It is useful to test the pseudoparticle method with our numerical 
example to assess its usefulness and accuracy. 

To introduce the pseudoparticle method, we rewrite Eq. (fT7|) in the form [3J 1 



f(rp,t) = [dr Q dp A(p-p Q + V r V(r,t)8t)8(r-r -^-8t)f(r ,p ,to), (46) 
J m 



where 



A(?r)= / (2W 6XP 



s-7T 2<ft ~ 1 /s-V n2 " +1 



(47) 



and 7T = p — p + V r V(r,t)5t. One notes that the second term in the exponential 
function is a sum involving third and higher powers of s. If one neglects these higher- 
order terms in the exponential function, then one obtains the approximation 

A(p -P + V r V(r, t)5t) » 6(p - p + V r V(r, t)8t) (48) 

and the time evolution of the Wigner function becomes [S] 

f(rp,t) « f L o(rp,t) = [dr dp 5(p-p + V r V(r,t)5t)5(r-ro-—8t)f{r ,p ,t$9) 

J m 

which leads to 

f(rp,t)af LO (rp,t)=f(r-p6t/m,p + V r V(r,t)6t, t ). (50) 

The above equation provides the basis for the lowest-order (LO) pseudoparticle 
approximation of the time evolution of the Wigner function One divides the phase 
space into cells (pseudoparticles) and follows the trajectories of these pseudoparticles 
using classical equations of motion, Eqs. (fT2j) and (JT3J). Eq. (J5U|) specifies that the 
Wigner function at the new phase space point (rp) at t — to + St is approximately the 
same as the initial Wigner function at the initial classical trajectory phase space point 
(r —p8t/m, p + V r V{r,t)5t) at to- The subscript LO in fLo( r P,t) is to indicate that it 
is the solution of the lowest-order pseudoparticle approximation obtained by following 
classical trajectories. This method was used in [Hj where it was found that the lowest- 
order pseudoparticle approximation reproduces quite well the general features of the 
time-dependent Hartree-Fock approximation. The use of classical trajectories to obtain 
an approximate time evolution of the Wigner function was also emphasized by Zachos 
and Curtright [T4] . 

How good is the lowest-order pseudoparticle approximation? The accuracy of the 
pseudoparticle approximation depends on the potential V(r,t). If the potential is a 
constant, a linear, or a harmonic oscillator potential, then even the LO pseudoparticle 
method gives the exact solution, as is clearly indicated in Eq. ()47j1 . The LO 
pseudoparticle solution deviates from the correct solution if the third and higher order 
spatial derivatives of the potential do not vanish. 
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Figure 2. Time dependence of the Wigner function obtained in different methods. 
The left panel gives the correct Wigner function obtained from eigenfunctions 
and eigenvalues, the middle panel gives the results from the lowest-order (LO) 
pseudoparticle approximation, and the right panel shows the Wigner function obtained 
by the next-to-lowest (NLO) order pseudoparticle approximation with ft 2 corrections. 



For the non-stationary state in a Gaussian potential in our test, we show the results 
of the LO pseudoparticle approximation in Figs. (2a'), (2b'), and (2c'), to be compared 
with the correct results from the eigenfunction method in Figs. (2a), (2b), and (2c). In 
our calculations, we specify the Wigner function at phase space coordinates of a fixed 
lattice. The Wigner function at the phase space points on the righthand side of Eq. 
(J5UJ) is obtained by interpolation using the Wigner function f(r ,p Q ,t ) at the earlier 
time to. As one observes, the general features of the oscillations are approximately 
reproduced. However, there are deviations from the correct results as time increases. 
For example, the Wigner function values for p — 0, t — 3, and \x\ ~ 2 obtained in the 
LO pseudoparticle approximation in Fig. (2a') are much larger than the corresponding 
correct values in Fig. (2a). The values of the Wigner function for p = 0.6, t = 3, and 
x ~ — 1 in the LO pseudoparticle approximation in Fig. (2c') are much smaller than the 
corresponding correct values in Fig. (2c). 

We can correct for the deviations of the LO pseudoparticle approximation. For a 
small time increment 5t, we can keep terms up to the first order in 5t in Eq. (|47jl . The 
function A (it) becomes 



A(7r) = / *L e i8 - n ' n 
1 2nh 



' 2zSt™ 1 fs-VV\ 2n+1 . ' 



(51) 
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The integration over s can be easily carried out and we obtain 



A(tt) 



2n 

1 + ^E(^TI)T(S] ^T +1 V(r,t) 



<J(tt), (52) 



where Vjf acts on V(r, t) and Vp acts on 5(ir). The corrected pseudoparticle solution 



7V 

. i I I I i I I i: i s < .-Jill I v 

P 

of the time evolution of the Wigner function is then 

2n 

\2n+l 



f(rp,t) 




V{r,t) 



ho(r,p,t), (53) 



where acts on fLo( r ,P,t)- It should be emphasized that no odd powers of h and 
even powers of (Vjf • V£) are present in the above expansion. We can improve upon 
the pseudoparticle approximation by including additional contributions in powers of h 2 
involving higher-order derivatives of V(r,t) and f'Lo{ r ,P,t)- 

We show in Figs. (2a"), (2b"), and (2c") the results in the next-to- lowest order 
approximation obtained by including the correction term of order h 2 , which involves 
the third derivatives of the potential and the Wigner function. In this calculation, we 
propagate the Wigner function from i to t = to + St for a small time increment St 
following classical trajectories as in Eq. ([49)1 . The Wigner function is then corrected 
at time t using Eq. ([53)1 . The corrected Wigner function at t is used in Eq. f!49|) to 
propagate to time t + St by following classical trajectories. It is then corrected at time 
t + St using Eq. These steps are repeated to obtain the time evolution of the 

Wigner function in this pseudoparticle approximation. To avoid numerical instability in 
the time evolution, we use the Wigner function fLo( r P,t) obtained in the lowest-order 
approximation to calculate the third momentum derivative, when we evaluate the h 2 
correction term in Eq. (|53|). 

Results in Fig. 2 indicate that the inclusion of the h 2 correction term in the next-to- 
lowest (NLO) order improves the pseudoparticle approximation. In particular, the NLO 
values of the Wigner function at p = 0, t = 3, and \x\ ~ 2 obtained with h 2 corrections 
Fig. (2a") are now closer to the values obtained with eigenf unctions. The NLO values 
of the Wigner function at p — 0.6, t = 3, and x ~ — 1 in Fig. (2c') are also closer to the 
corresponding values obtained with eigenfunctions. Although small deviations remain 
in the NLO solution, the general features of the Wigner function are reasonably well 
reproduced. 

We give here an alternative method to evaluate the correction terms in the 
pseudoparticle method. As one notes in Eq. (|52|). the correction terms involve high- 
order derivatives of the 5-function. They can be evaluated approximately by expressing 
the 5-function in terms of known functions. Making use of the othonormality of the 
harmonic oscillator wave functions {</> m (a;)}, we can represent a one-dimensional delta 
function S(x) by 

00 2 2 00 ( —'\\ m 

= £ C(*)<U0) = ~^e-^ £ H 2m {ax))—l- (54) 

m=0 V 71 " m=0 Z m - 

where H2 m (ax) is the Hermite polynomial of order 2m, and a is a width parameter. In 
practice, the summation in the above equation is carried out up to an order M, and the 
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delta function 5(x) can be approximated by the distribution 

5{x) « D(x) = A,®e-^ £ (55) 
v 71 " m=0 2 m - 

where Am is 

M (2m)! 



.4 



M 



m _ n 2 2m m!m! 



(56) 



and t4a/ is so chosen that D(x) is normalized to J D(x)dx = 1. A three-dimensional 
generalization of this D-function, D(x), can be obtained as the product of three one- 
dimensional /^-functions of their respective component coordinates, D(x x )D(x y )D(x z ). 
Using a distribution function of this type, we can write the explicit solution of the 
Wigner function as 

f{rp,t) = ho(r,p,t) + St fdr dp S(r - r - —6t)f(r ,p Q ,t ) 

J m 



2n 

X E (2^Tl)! (jU {(y?-V?f n+1 V(r,t)D(p-p + StV r V(r,t))\, (57) 

where applies to D{jp — p + 5tV r V(r, t)). As D is a known analytical function 
of p, (Vp) 2n+1 D(p — p + StV r V(r,t)) can be obtained analytically. The higher-order 
derivatives of the Wigner function in Eq. can be converted into an integration over 

Po- 

The pseudoparticle method utilizes classical trajectories and is easy to use. It is 
even exact for a constant, linear, or harmonic oscillator potential. The development of 
systematic higher-order corrections to the pseudoparticle method for a general potential 
in Eq. (f53|) allows one to improve on the approximate solution to make the pseudoparticle 
method a useful tool for future applications. 

In applying the pseudoparticle method, one can use the Eulerian picture with a 
fixed lattice or alternatively the Lagrangian picture with a set of pseudoparticle phase 
space coordinates, as in molecular dynamics. The distribution D of Eq. (J55|) can be 
used to go from one picture to the other. For example, if one starts with a Wigner 
function in the Eulerian picture, one divides the phase space into pseudoparticle cells 
with volume element AviAp { and specifies the phase space coordinates (r*^) and its 
Wigner function /(r^pj. The set of variables (r^) and /i(T*iP i )(= /(r^pj) can be 
used in the Lagrangian picture to evolve the phase space dynamics. Conversely, if one is 
given the set of variables (r^pj and /r(j'iPj) in the Lagrangian picture, then the Wigner 
function in the Eulerian picture /e(t*p) at coordinates (rp) is given by 

f E (rp) = ]T D r (r - r t )D p (p - Pi )h{r iPi ) Ar.APi, (58) 

i 

where the width parameter a in D r (and D p ) depends on the magnitude of Ar^ (and 
Apj). These transcriptions between the two pictures facilitate the application of the 
pseudoparticle method. 



Explicit Solution of the Time Evolution of the Wigner Function^ 



15 



6. Conclusions and Discussions 

We have reviewed and applied here the explicit solution of the time evolution of the 
Wigner function obtained previously in 1982 ||. The basic idea is to represent the 
Wigner function in terms of auxiliary phase space coordinates, which obey simple 
equations of motion. These equations are similar to the classical equations of motion. 
They can be solved easily. The solutions of these equations of motion can then be 
used to evaluate the time evolution of the Wigner function. We have demonstrated the 
usefulness of the explicit solution using a numerical example. We find that the explicit 
solution leads to the correct time evolution of the Wigner function, even for a Wigner 
function with strong spatial and temporal variations and regions of negative values. 

We have also reviewed and tested the pseudoparticle method for the evaluation of 
the time evolution of the Wigner function. For our example of a non-stationary state in 
a Gaussian potential, the lowest-order pseudoparticle approximation gives the correct 
features of the time evolution, but there are deviations from the correct results. We have 
developed a systematic way to improve the pseudoparticle method involving correction 
terms in powers of h 2 containing high-order derivatives of the potential and the Wigner 
function. 

The simplicity of the different methods discussed here will facilitate their 
application to quantum dynamics in phase space. They can be used to study quantum 
particle dynamics in a time-dependent, multi-dimensional potential. They can also be 
applied to study the one-body Wigner function of a many-particle system in a time- 
dependent mean-field potential, as in With the addition of a collision term, they 
can be used to describe the dynamics of a quantum Boltzmann equation. The explicit 
solution of Eqs. (fTT)|) - (|2"Tf or (pSjl - P^j) is probably best handled in a lattice of phase 
space points, as the sine and cosine transforms are simplest in such a lattice. On 
the other hand, the pseudoparticle method of Eqs. and can be investigated 
both in terms of pseudoparticle phase space coordinates in the Lagrangian picture or 
alternatively in a fixed lattice of phase space points in the Eulerian picture. The use of 
pseudoparticle coordinates may provide substantial saving of computer storage capacity 
when a significant fraction of the phase space is empty. Future research using these 
methods will allow us to explore further the richness of quantum dynamics in phase 
space, which Professor Wigner first pioneered for us. 
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